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Abstract 

The present work proposes a theory of isotropic and homogeneous turbulence for incompressible 
fluids, which assumes that the turbulence is due to the bifurcations associated to the velocity field. 
The theory is formulated using a representation of the fluid motion which is more general than 
the classical Navier-Stokes equations, where the fluid state variables are expressed in terms of the 
referential coordinates. 

The theory is developed according to the following four items: 1) Study of the route toward the 
turbulence through the bifurcations analysis of the kinematic equations. 2) Referential description 
of the motion and calculation of the velocity fluctuation using the Lyapunov analysis of the local 
deformation. 3) Study of the mechanism of the energy cascade from large to small scales through the 
Lyapunov analysis of the relative kinematics equations of motion. 4) Determination of the statistics 
of the velocity difference with the Fourier analysis. Each item contributes to the formulation of 
the theory. 

The theory gives the connection between number of bifurcations, scales and Reynolds number at 
the onset of the turbulence and supplies an explanation for the mechanism of the energy cascade 
which leads to the closure of the von Karman-Howarth equation. The theory also gives the statistics 
of the velocity difference fluctuation and permits the calculation of its PDF. 

The presented results show that the proposed theory describes quite well the properties of the 
isotropic turbulence. 

PACS numbers: Valid PACS appear here 
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I. INTRODUCTION 



This work presents a theory of isotropic and homogeneous turbulence for an incompress- 
ible fluid formulated for an infinite fluid domain. The theory is mainly motivated by the 
fact that in turbulence the fluid kinematics is subjected to bifurcations Ij and exhibits a 
chaotic behavior and huge mixing [2I, resulting to be much more rapid than the fluid state 
variables. This characteristics implies that the accepted kinematical hypothesis for deriving 
the Navier-Stokes equations could require the consideration of very small length scales and 
times for describing the fluid motion [s^ and therefore a very large number of degrees of 
freedom. To avoid the difficulties arising from the consideration of these small scales, the 
referential description of motion is adopted, where the fluid state variables are expressed in 
terms of the so called referential coordinates which coincide with the material coordinates 
for a given fluid configuration 3|. 

The other very important subjects of the turbulence are the non-gaussian statistics of the 
velocity difference and the mechanism of the kinetic energy cascade. This latter is directly 
related to the relative motion of a pair of fluid particles jjO, I^, Q] and is responsible for the 
shape of the developed energy spectrum. 

For these reasons the present theory is based on: 



1. Landau hypothesis, following which the turbulence is caused by the bifurcations of the 
velocity field 

2. Referential description of motion, where velocity field and stress tensor are mapped 
with respect to the referential coordinates j^. 

3. Study of the energy cascade through Lyapunov analysis of the relative kinematics. 

4. Statistical analysis of the velocity difference fluctuations. 

In the first part of the work, the road toward the turbulence is studied through the 
bifurcations analysis of the kinematic equations. These bifurcations arise from the math- 
ematical structure of the velocity field, where the Reynolds number plays the role of the 
"control parameter". This analysis supplies the connection between number of bifurcations 
and the critical Reynolds number for isotropic turbulence, showing that the length scales 
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are continuously distributed and that each of them is important for the description of the 
motion. 

In the second part, the momentum equations are formulated according to the referential 
representation of motion, whereas the kinematics of the local deformation is studied with 
the Lyapunov theory. The fluid motion is described adopting the referential configuration 
which corresponds to the fluid placement at the onset of this fluctuation. This choice allows 
the velocity fluctuations to be analytically expressed through the Lyapunov analysis of the 
kinematics of the fluid deformation. 

The third part deals with the relative kinematic between two trajectories, which is also 
analyzed with the Lyapunov theory. This analysis gives an explanation of the mechanism 
of kinetic energy transfer between length scales and leads to the closure of the von Karman- 
Howarth equation [g] (see Appendix), where the unknown function K{r), which represents 
the inertia forces, is here expressed in terms of the longitudinal correlation function. The 
obtained expression of K{r) satisfies the conservation law which states that the inertia forces 



only transfer the kinetic energy 

To complete the theory, the statistics of velocity difference is studied through the Fourier 
analysis of the velocity fluctuations. An analytical expression for the velocity difference 
and for its PDF is obtained in case of isotropic turbulence. This expression incorporates 
an unknown function, related to the skewness, which is immediately identified through the 
obtained expression of K{r). 

Finally, the several results obtained with this theory are compared with the data exist- 
ing in the literature, indicating that the proposed theory adequately describes the various 
properties of the turbulence. 

II. BIFURCATION ANALYSIS OF THE KINEMATIC EQUATIONS 

In this session, the route toward the turbulence is studied through the analysis of the bi- 
furcations of the kinematic equations. To analyze this question, a viscous and incompressible 
fluid in the infinite domain is considered, whose kinematic equations are 

^ = u(x,t;i?e) (1) 
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where x and Re are the position and Reynolds number, whereas u(x, t; Re) is a single 
realization of the ensemble of the velocity fields, written in the reference frame 3?, which 
satisfies the Navier-Stokes equations 



V-u = 



— + uVu + — - t/V^u = 
ot p 



(2) 



p and u are, respectively, density and kinematic viscosity whereas p is the fluid pressure 
which can be eliminated by taking the divergence of the momentum equation ^ 



+ Vu : Vu = 



(3) 



Now, let consider an assigned velocity fleld at a given time, and the flxed points X of Eq. 
([T]) which satisfy to rfX/rft = 0. Increasing the Reynolds number, X will vary according to 
Eq. ([T]), which can be solved by the continuation method j^, 



X = Xn 



Vu^^TT^ dRe 

Reo 



dRe 



(4) 



where Xq is the flxed point calculated at Re = Re^. The Reynolds number influences the 
mathematical structure of Eq. ([T]) through the Navier-Stokes equations in such a way that, 
for small i?e, the viscosity forces which are stronger than the inertia ones, make u an almost 




Re^ Rep 



FIG. 1: Map of the bifurcations. 
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smooth function of X. When the Reynolds number increases, as long as the Jacobian Vu is 
nonsingular, X exhibits smooth variations with i?e, whereas at a certain i?e, this Jacobian 
becomes singular due to the higher inertia- viscous forces ratio, resulting det (Vu) = 0. This 
can correspond to the first bifurcation, where at least one of the eigenvalues of Vu crosses 
the imaginary axis and X appears to be discontinuous with respect to Re j^, . Increasing 
again the Reynolds number, X will show smooth variations until to the next bifurcation. 

Figure 1 shows a scheme of bifurcations, where the component X of X is reported in 
terms of Reynolds number. Starting from Rcq, the diagram is regular, until to Rep, where 
the first bifurcation determines two branches, whose distance AXp is measured at the next 
bifurcation. For each bifurcation, AX gives a length scale of the velocity field at the current 
Reynolds number, whereas ARe represents the distance between two successive bifurcations. 
After P, Eq. (jlj) does not indicate which of the two possible branches the system will choose, 
thus a bifurcation causes a lost of informations with respect to the initial data [10] . Therefore, 
the fiuctuations are important for the choice of the branch that the system will follow 

Further increments of Re cause an increment of the number of bifurcations whose scaling 

^ _ ^ _ A Ren 

ARen+1 

For Re oo, the convergence of a„ and (5„ is not granted in general, whereas for period- 



l3- 



doubling bifurcations, these admit the following limits H] 



a = lim |a„| = 2.502... 5 = lim |(5„| = 4.669... (6) 

Re—>oo Re—foa 

These are the famous Feigenbaum numbers, which are two universal constants, independent 
on the mathematical details of the period-doubling bifurcations. For bifurcations of other 
kind, an and 6n can converge to different values or can oscillate around to average values. 

In the present analysis, the length scales In = AXn are assumed to be expressed by the 
asymptotic approximation 

ln = -K (7) 

Equation (^^) supplies the length scales in terms of the numbers of the bifurcations encoun- 
tered along a given path of fixed points, where a is the Feigenbaum constant given by Eq. ([6]) 



and li represents the maximum length scale. According to 



1^ . the bifurcations 



generate a route toward the chaos which depends on n. As long as n < 2, each bifurcation 
adds a new frequency into the power spectrum of u and this corresponds to hmit cycles or 
quasi periodic motions, whereas for n > 3, the situation drastically changes, since u exhibits 
more numerous frequencies and this generates chaotic motion 13] • This occurs for a 
single realization of the ensemble of the velocity field. The fluctuations of u(x, t) will cause 
further variations of the several scales /„ in Eq. ([7]), thus the bifurcations maps will be more 
complicated than Fig. [H and the recognizing the diverse scales and bifurcations could not 
be possible. This is a scenario with continuously distributed length scales, where all of them 
are important for describing the fluid motion. 



A. Critical Reynolds number 

Equation (JTj) describes the route toward the chaos and is assumed to be valid until the 
onset of the turbulence. In this situation the minimum for /„ can not be less than the 
dissipation length or Kolmogorov scale £ = [v'^ /sY^'^ [l], where e is the energy dissipation 
rate (see Appendix), whereas li gives a good estimation of the correlation length of the 
phenomenon j^, Q which, in this case is the Taylor scale At- Thus, £ < In < At, and 

" (8) 



where is the number of bifurcations at the beginning of the turbulence. 

Equation ([8]) gives the connection between the critical Reynolds number and number of 
bifurcations. In fact, the characteristic Reynolds numbers associated to the scales i and At 
are Rk = iux/i^ = 1 and Rx = Xtu/u, respectively, where uk = {veY^^ is characteristic 
velocity at the Kolmogorov scale, and u = a/ (uiUi) /3 is the velocity standard deviation [3]. 
For isotropic turbulence, these scales are linked each other by 



a 



At/^=15^/^/rI (9) 
In view of Eq.®, this ratio can be also expressed through A^. 



a 



(10) 



The value Rx ~ 1.613 obtained for A^ = 2 is not compatible with At which is the correlation 
scale, while the result Rx — 10.12, calculated for A^ = 3, is an acceptable minimum value 



"or R\. This result agrees with the various scenarios describing the roads to the turbulence 
~ 3, and with dive.e expen^ent. UUU wMch .a. tha. turbulence 



Ol - TLA- J-i 



begins for > 3. Of course, this minimum value for Rx is the result of the assumptions 



a 



2.502, h ^ At, / 



N 



and of the asymptotic approximation ([7]). 



III. REFERENTIAL DESCRIPTION OF MOTION. VELOCITY FLUCTUATION 



Now, we present a formulation of the fluid equations of motion which is based on the 
referential description of the motion. This formulation is more general than the classical 
Navier-Stokes equations and is capable to take into account the effects of the fluid kinematics 
which can be much faster than the fluid state variables. This description of motion allows to 
calculate the velocity fluctuation through the Lyapunov analysis of the local deformation. 

This representation of motion is based on the fact that a given fluid property Q is an 
explicit function of the referential displacement xq and of the time j^, i.e. 

n = n{^o,t) (11) 

The referential displacement coincides with the material position for a given fluid configura- 
tion, thus xo plays the role of the label which identifies the specific fluid particle Since 

Current Configuration 
(Deformed) 





Reference Configuration , ' 
(Undeformed) 

t=o S^, 



FIG. 2: Referential description of the kinematics of deformation. 
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any fluid motion has infinitely many different referential descriptions which are equally valid 
[sl, it is convenient to choose the referential configuration corresponding to the fluid place- 
ment at the onset of the deformation (see Fig. [2]). According to Truesdell js], n(xo,t) and 
its derivatives with respect to xq are supposed to be smooth functions of t and xq. Hence, if 
X = x(xo, t) represents the fluid motion, Q is expressed in terms of the geometrical position 
X, through the inverse of x, xq = x^^(x, t) 

fi(x,t) = fi(x-'(x,t),t) (12) 



and its derivative with respect to x is 



(13) 



9x 9xo 9x 

The bifurcations of Eq. ([1]) make x ^ singular transformation, thus, in proximity of a 
bifurcation, dVL/d^ varies much more quickly than dVt/d^Q because of the local stretching 
9x/9xo, which is here calculated with the Lyapunov theory, as 

1^ ^ e^* (14) 

9X0 

where A = max(Ai, A2, A3) is the maximal Lyapunov exponent and Aj, (i = 1, 2, 3) are the 
Lyapunov exponents. Due to the incompressibility, Ai + A2 + A3 = 0, thus, A > 0. 

The velocity fluctuation of the particle xq -or Lagrangian fluctuation- is calculated using 
the momentum equations, where stress tensor and velocity field are mapped with respect to 
the referential coordinates at the beginning of the deformation 

duk\ ldTkh,__ ,^dxjo 



where (duk/dt)^^ is the acceleration of the particle xg, whereas T^h represents the stress 
tensor 

^ . f dufi duk\ ,^ . 

Thk = -pShk + U— + ^ (16) 



dxk dxh, 

Note that, Eq. f|T5|) is more general than the classical Navier-Stokes equations, since it can be 
applied to fluid particles which exhibit non-smooth displacements and irregular boundaries 
31], as in the present case. Since 9x/9xo is much more rapid than dTkh/dxjo, this fluctuation 
is calculated integrating Eq. (ITSl) from t = to 00, considering dTkh/dxjo constant with 
respect to (9x/c?xo, i.e. 
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The velocity fluctuation in a fixed point of space x -or Eulerian fluctuation- is calculated 
taking into account the expression of the Eulerian time derivative of Uk, which is js] 



duk 
dt 



duk 
dt 



duk , ^.dxho 

^0,t)— Ui 



Therefore, this velocity fluctuation is 

UkU) 



1 

A 



dx 



hO 



dxj 



duk 
dt 



(18) 



(19) 



These velocity fluctuations, which stem from the bifurcations of the velocity field, do not 
modify the average values of the momentum and of the kinetic energy of fluid. 



IV. LYAPUNOV ANALYSIS OF THE RELATIVE KINEMATICS 

In order to investigate the mechanism of the energy cascade, the properties of the relative 
kinematic equations are here studied with the Lyapunov analysis. These equations are 

^ = u(x,t), ^ = u(x',t) (20) 

where u(x, t) = (ui,M2,M3), u(x',t) = u' = {u'l.u^.u'^, whereas ui and u[ are the velocity 
components expressed in the reference frame 3?. Since the bifurcations do not modify the 
total momentum and kinetic energy, the solutions of Eq. fl20|) preserve these quantities. 
With reference to Fig. [31 these solutions correspond to the paths, x(t) and x'(t), located 



into a material volume ^(t) which changes its geometry according to the fluid motion |18| . 
whereas its volume remains unaltered. This is a toroidal volume, where 5*^ and R are, 
respectively, the poloidal surface and the toroidal dimension of E, whereas X and X' are the 
intersections of x(t) and x'(t) with 5*^, where r = |X' — X| is the poloidal dimension, thus 
Sp ^ r"^. The velocity difference components Am„ = u'^ — Un and Am^ = u'^ — Ur lay on Sp 
and are normal and parallel to r, respectively, whereas Ub is the average velocity component 
along the direction normal to Sp. The equations describing the evolution of these quantities 
preserve the volume and the momentum of S. These can be written as 

Jt ^^'^^ = ° ^^^^ 
Jt = ° 

(22) 

d , 



FIG. 3: Scheme of the relative kinematics of two fluid particles 



Equations (12T]) and (l22l) represent, respectively, the continuity equation, and the momentum 



18|- 



equations according to the third Helmholtz theorem on the vorticity 

The Lyapunov analysis, applied to Eqs. ([21]) and (122]) . states that R ^ Rq e^*, hence, 
Eqs. (I2T]) and ^ become 



^ = XAul (23) 

dub , 
^ = -^"^ 

where A(r) >0 is the maximal finite scale Lyapunov exponents associated to Eqs. (120]) . with 
A(0) = A. As the result, -ut and Au ^ Am„ oc e'^/^*. 
Now, it is worth to remark that the following quantity 

T.l(u.u') (24) 
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expresses the transfer of the kinetic energy between the points x and x'. Its average (T) is 
calculated on the ensemble of the diverse pairs of trajectories which pass through X and X' 
and which are contained into the various toroidal volumes. This average is obtained from 
Eqs. fl2^ . taking into account the homogeneity, the isotropy and the time independence 
upon the time of the average kinetic energy ( d{u ■ u.)/dt = 0). 

= <1 E = ^^'i9-f) (25) 

i=r,n,b 

f and g are longitudinal and lateral velocity correlation functions, that, because of the 
incompressibility, are related each other through Eq. fIBUl) (see Appendix). Thus, (T) is 

(TT) = Kr)r (26) 

If T were an ergodic function, its average on the statistical ensemble should coincide with 
the average over time which in turn is equal to zero since T is the time derivative of u ■ u'. 
As the consequence, there would not be any transfer of energy between the parts of fluid. 
Therefore, the fluid incompressibility is a sufficient condition to state that T is a non ergodic 
function, whose statistical average is determined as soon as A is known. To calculate A, it is 
convenient to express the velocity difference Au = u(x', t) — u(x, t) in the Lyapunov basis 
E associated to Eqs. (l20l) . which is made by orthonormal vectors arising from Eqs. (l20l) 
jiol . Q. The velocity difference expressed in E, Av = {v[ — vi, V2 — f2, v'^ — V3), satisfies 
the following equations, which hold for t — > cxd 

w- -Vi = Xi fi, z = 1, 2, 3 (27) 

where fj, Vi and v'^ are, respectively, the components of r = x' — x, u(x, t) and u(x',t) 
written in E. Then, Aur and r can be expressed in terms of Av and f as 

r ^ ^ ■ Qr, Aur^^- QAv (28) 

Into Eqs. fl28l) . Q is the fluctuating rotation matrix transformation from E to 3?, and 
^ = (X' — X)/|X' — X|. The standard deviation of Aur is calculated from Eqs. (l28l) . taking 
into account the isotropy and that Av ^ Af 

{Aulir)) = AV^ (29) 
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This standard deviation can be also expressed through the longitudinal correlation function 
/ 

{Aulir)) = 2u\l-fir)) (30) 

being u the standard deviation of the longitudinal velocity. The maximal Lyapunov exponent 
is calculated in function of /, from Eqs. fl29|) and fl30l) 



A(r) = -v/2(l-/(r)) (31) 



Hence, substituting Eq. (131]) into Eq. (!26l) . one obtains the expression of (T) in terms of 
the longitudinal correlation function 



2 dr 

where, thanks to the isotropy, (T) is a function of r alone 



m=A/VI^ (32) 



Equation 



reflects the well known property of the inertia forces of transferring the 



kinetic energy 7(1 between the several regions of the fluid domain. 

V. CLOSURE OF THE VON KARMAN-HOWARTH EQUATION 

The closure of the von Karman-Howarth equation is now carried out using the previous 
Lyapunov analysis. 

The function K{r) is defined through the following relation (see also Eq. flB2|) in the Ap- 
pendix) 

Mu, - 4)) = + 3K{r) (33) 

The repeated indexes into Eq. ( l33l) . i and k, indicate the summations with respect to the 
same indexes. In order to obtain the expression of K{r), it is worth to remark the following 
identity 

Uiu'iiuk - u'f^) = T rj, - -^{uiu'-rk) (34) 

The average of Eq. (IMI) is calculated on the ensemble of the trajectories passing through 
X and X'. It is supposed that the ergodic hypothesis holds for the last term at the right 
hand-side of Eq. (!34l) . thus this latter can be calculated through the average over time. Since 
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this term is the time derivative of Uiu'^Vk, this gives null contribution. Hence, accounting for 
the isotropy, one obtains 

J-(n.n:K-n',)) = ^r + 3(T) (35) 



Comparing Eqs. (155]) and and taking into account that K{0) = K{r) = (T), i.e. 



Equation fl36|) represents the proposed closure of the von Karman-Howarth equation, and 
expresses the transfer of kinetic energy between the diverse fluid regions. This is a kinematic 
mechanism, caused by the bifurcations cascades of Eq. (l20i) . which preserves total momen- 
tum and kinetic energy. The analytical structure of Eq.( l36l) states that this mechanism 
consists of a flow of the kinetic energy from large to small scales which only redistributes 
the kinetic energy between wavelengths. 



The skewness of Aur is determined once K{r) is known This is 

The longitudinal triple correlation k{r) is calculated by Eq. (163!) (see Appendix). Since / 
and k are, respectively, even and odd functions of r with /(O) = 1, A;(0) = k'{0) = k"{0) =0, 
H^iO) is given by 

^3(0) = lim/73(r) = p^p^ (38) 

where the apex denote the derivative with respect to r. To obtain H^lO), observe that, near 
the origin, K behaves as 

K = «3yr7^/"(o)^ + O(r^) (39) 

then, substituting Eq. fl39l) into Eq. (|63|1 (see Appendix) and accounting for Eq. fl38l) . one 
obtains 

H3{0) = -^ = -0.42857... (40) 
This value of H^lO) is a constant of the present theory, which does not depend on the 



Reynolds number 



nn 



his is in agreement with the several sources of data existing in the 



literature such as [7|, l2l|, |22|, |23|| (and Refs. therein) and the knowledge of it gives the entity 



of the mechanism of energy cascade. 
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VI. STATISTICAL ANALYSIS OF VELOCITY DIFFERENCE 



Although the previous analysis leads to the closure of the von Karman-Howarth equation, 
it does not give any information about the statistics of velocity difference Au(r) = u(x + 
r) - u(x). 

In this section, the statistical properties of Au(r), are investigated through the Fourier 
analysis of the velocity fluctuation given by Eq. f|T9l) . This fluctuation is 

u = $:U(,.)e---i$:^(..)e-- (41) 
where U(k) = {Ui{k,), U2{k), [/^{k)) are the components of velocity spectrum, which satisfy 

3 

(42) 

All the components U(k;) ~ d\J{K,)/dt/A are random variables distributed according to 
certain distribution functions, which are statistically orthogonal each other 7|. 

Thanks to the local isotropy, u is sum of several dependent random variables which are 
identically distributed therefore u tends to a gaussian variable 24|, and U(k;) satisfies 
the Lindeberg condition, a very general necessary and sufficient condition for satisfying 
the central limit theorem [2^. This condition does not apply to the Fourier coefficients of 
Au. In fact, since Au is the difference between two dependent gaussian variables, its PDF 
could be a non gaussian distribution function. In x = 0, the velocity difference Au(r) = 
(Ami, Au2, AM3) is given by 

A«,^l E ^(e^'^- -l)^L + 5 + P + iV (43) 



This fluctuation consists of the contributions appearing into Eq. fl42|) : in particular, L 
represents the sum of all linear terms due to the viscosity and B is the sum of all bilinear 
terms arising from inertia and pressure forces. P and are, respectively, the sums of 
definite positive and negative square terms, which derive from inertia and pressure forces. 
The quantity L + B tends to a gaussian random variable being the sum of statistically 
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orthogonal terms 2J, |25|, while P and do not, as they are linear combinations of squares 



251]. Their general expressions are (25 1 

P = Po + ?7i + Vl 

(44) 

AT = + Ci + a 

where Pq and A^^q are constants, and rji, ri2, Ci aiid ^2 are four different centered random 
gaussian variables. Therefore, the fluctuation Aup with zero average reads as 

An, = Mr)^ + Mr) {xiv' " 1) " (C' " 1)) (45) 

where ^, rj and ( are independent centered random variables which have gaussian distribution 
functions with standard deviation equal to the unity. The parameter x is a function of 
Reynolds number, whereas ipi and ■?/'2 are functions of space coordinates, which also depend 
on the Reynolds number. 

At the Kolmogorov scale the order of magnitude of the velocity fluctuations is ux'^r/i, 
with r = 1/A, and tp2 is negligible because is due to the inertia forces: this immediately 
identifies ipi ~ uk^t/ I. 

On the contrary, at the Taylor scale, ipi is negligible and the order of magnitude of the 
velocity fluctuations is v?t/\t-, therefore 'ip2 ~ u^t/Xt and the ratio is a function of 

Rx 



V,(r,fl,) = ^«4^=. (46) 

^IJiir) uk^Xt y 15V15 

where ip{r,Rx) = 0(1), is a function which has to be determined. Hence, the longitudinal 
velocity difference Aur, is written as 

A^. ^ e + ^(x(r/^-l)-(C^-l)) 

The quadratic term at the right hand side of Eq. ( H71) represents the velocity fluctuations at 
the bigger scales, and there is no physical reason for which this must be bounded between 
same limits. Consequentely, x must be a definite positive function of Pa- 

Equation fH7|) gives the mathematical structure of Aur, whose dimensionless statistical 
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moments are easily calculated considering that ^, t] and ( are independent gaussian variables 

^ _ (A<) 1 



EQ^'r~')((x(^'-i)-(c^-i))') 



(48) 



where 



(49) 



1=0 ^ ^ 

In particular, the third moment or skewness, H^, which is responsible for the energy cascade, 
is 

H,^^tit^^ (50) 
(1 + 2'*' (1 + X-')P 

For X 7^ 1) the skewness and all the odd order moments are different from zero, and for n > 3, 
all the absolute moments are rising functions of Rx, thus Am,, exhibits an intermittency 
whose entity increases with the Reynolds number. If and x were both known, the other 
statistical moments can be consequentely calculated with Eq. (HHj) . The function ip{r,Rx) 
is determined for Aur from Eqs. ( l50ll and ( !37l) . For r=0, one obtains the relationship 

8^0^ (1 - X') ^ 3 
(l + 2^o^(l + X^))^/^ 7 

ipo = ip{^,Rx) =0(1); is given by Eq. (1461) . where, its exact value has to be calculated, 
whereas x is a positive function of Rx which must also be determined. To determine such 
quantities, note that Eq. flSlI) is an algebraic relationship which gives x terms of Rx, as 
shown in Fig. HI In any case, x exhibits the limit x — 0.86592 for Rx oo, whereas Rx 
admits the minimim {Rx)min which depends on "ipQ. Below such minimum, Eq. f l5T|) does 
not admit solutions with x > 0. Then, according to the analysis of section Hi Al ipQ is chosen 
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FIG. 4: Parameter x plotted as the function of Rx. 



in such a way that (-Ra) 



mm 



10.12 as shown in Fig. |U resulting i/jq ~ 1.075. Now, all the 



moments of Aur can be calculated by Eqs. (HH!) and (149|) in terms of -Ra- 
The PDF of Aur is expressed through the Frobenious- Perron equation 

F{Au'r) = [ [ [piOpiv)piO S (Aur-Au'r) d^drjdC (52) 

where Aur is calculated with Eq. ( 1471) . S is the Dirac delta and p is a gaussian PDF whose 
average value and standard deviation are equal to and 1, respectively. 

For non-isotropic turbulence or in more complex cases with boundary conditions, the 
velocity spectrum could not satisfy the Lindeberg condition, thus the velocity will be not 
distrubuted following a Gaussian PDF, and Eq. (H5l) changes its analytical form and can in- 
corporate more intermittant terms which give the deviation with respect to the isotropic 
turbulence. Hence, the absolute statistical moments of Aur will be greater than those cal- 
culated with Eq. P7!) . indicating that, in a more complex situation than the isotropic 
turbulence, the intermittancy of Aur can be significantly stronger. 

VII. RESULTS AND DISCUSSION 

The results calculated with the proposed theory are now presented. 
As the first result, the evolution in time of the correlation function is calculated with 
the proposed closure of the von Karman-Howarth equation (Eq. (I36l) ). where the boundary 
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conditions are given by Eq. (165|) . The turbulent kinetic energy and the spectrums E{k,) and 
T(k) are calculated with Eq. fl66l) and Eqs. fl67|) . respectively. The calculation is carried 
out for the initial Reynolds number of Re = u{0)Lr/i' = 2000, where and u{0) are, 
respectively, the characteristic dimension of the problem and the initial velocity standard 
deviation. The initial condition is a gaussian correlation function with Xt/L^ = 1/(2^/2). 
The dimensionless time of the problem is defined as i = t u{0)/Lr. 

Equation fl6T]) was numerically solved adopting the Crank-Nicholson integrator scheme 
with variable time step, where the discretization of the space domain is made by — 1 
intervals of the same amplitude Ar. This corresponds to a discretization of the Fourier 
space made by — 1 subsets in the interval [0, km], where = tt/ (2Ar). For the adopted 
initial Reynolds number, the choice N = 1500, gives an adequate discretization, which 
provides Ar < i, for the whole simulation. During the simulation, T(/t) must identically 
satisfy Eq. (!68!) (see Appendix) which states that T(k) does not modify the kinetic energy. 
To verify Eq.( !68|) . the integral of T{k,) is calculated with the trapezes rule from until to 
K,M, at each time step, therefore, the simulation will be considered to be accurate as long as 

poo 

/ T{K)dK ^ / T{K)dK = (53) 
Jo Jo 

namely, when the energy is distributed for k < km- As the simulation advances, according 
to Eq. (1361) . the energy cascade determines variations of E{k) and T(/t) at the higher wave- 
numbers, then Eq. ( l53l) will hold until to a certain time. For this reason, the simulation is 



stopped as soon as the following condition is achieved 
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T{K)dK\ > / \T{K)\di^ (54) 



N^Jo 

At the end of several simulations, we have Ar ^ 0.8 i, and, in this situation, the energy 
spectrum is here considered to be fully developed. 

The diagrams of Fig. [5]show the correlation functions /(r) and k{r) vs. the dimensionless 
distance r/A^, at different times of simulation. The kinetic energy and Taylor scale diminish 
according to Eqs. (|36l) and (1661) . thus /(r) and k{r) change in such a way that the length 
scales associated to their variations diminish as the time increases, whereas the maximum 
of |fc| decreases. At the final instants of the simulation, one obtains that / — 1 = 0( r^/^) 
for t/Xt = 0(1), whereas the maximum of \k\ is about 0.05. These results are in very 
good agreement with the numerous data of the literature I?! which concern the evolution of 
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FIG. 5: Correlation functions, / and k versus the separation distance at the times of simulation t 
= 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.63. 

correlation function and energy spectrum. Figure [6] shows the diagrams of E{k) and T{k) for 
the same times, where the dashed line in the plot of E{k), represents the —5/3 Kolmogorov 
law j^. The spectrums E{k) and T{k) vary according to Eqs. (!36l) and (!67|) . and, at the 
end of simulation, E{k) is about parallel to the dashed line in an opportune interval of the 
wave-numbers which defines the so called inertial range of Kolmogorov. This arises from the 
developed correlation function, which behaves like / — 1 = O (r^/^) for r = 0{Xt)- 

Next, the Kolmogorov function Q{r) and Kolmogorov constant C, are determined with 
the proposed theory, using the previous results of the simulation. 

Following the Kolmogorov theory, the Kolmogorov function, which is defined as 



Q{r) 



re 



(55) 



is constant with respect to r, and is equal to 4/5 as long as r/Ay = 0(1). As shown in Fig. 
[71 for i = 0, the maximum of Q{r) is much greater than 4/5 and its variations with t/At 
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FIG. 6: Plot of E{k) and T{k) at the diverse times of simulation. 
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FIG. 8: (a) Maximum finite size Lyapunov exponent at the times of simulation t = 0, 0.1, 0.2, 0.3, 
0.4, 0.5, 0.6, 0.63; (b) and (c) skewness and Flatness versus r/Ar at t = and t = 0.6, respectively. 

can not be neglected. This is due to the arbitrary choice of the initial correlation function. 
At the successive times, the variations of / determine that the maximum of Q{r) and its 
variations decrease until to the final instants, where, with the exception of r/Ar ~ 0, Q{r) 
exhibits a qualitatively flat shape in a wide range of t/Xt, with a maximum which is quite 
close to 0.8. 

The Kolmogorov constant C is also calculated by definition 

2/3 



E(k) = C- 



K 



5/3 



(56) 



This is here determined, as the value of C which makes the curve represented by Eq. (156]) 
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FIG. 9: PDF of the velocity difference fluctuations at the times t =0 (a), t = 0.5 (b) and t =0.6 
(c). Continuous lines are for r =0, dashed lines are for r/Xx =1, dot-dashed lines are for r/Ar =5, 
dotted lines are for gaussian PDF. 

to be tangent to the energy spectrum E{k) previously calculated. At end simulation, C ~ 
1.932, namely C and Qmax agree very well to the corresponding quantities known from the 
literature. For the same simulation. Fig. [8^ shows the maximal finite scale Lyapunov expo- 
nent, calculated with Eq. (13T1) . where A varies according to /. For t = 0, the variations of A 
are relatively small because of the adopted initial correlation function which is a gaussian, 
whereas as the time increases, the variations of / determine sizable increments of A and of 
its slope in proximity of the origin. Then, for developed spectrum, since / — 1 = 0(r^/'^), the 
maximal finite scale Lyapunov exponent behaves like A ~ Thus, the diffusivity coef- 

ficient associated to the relative motion between two fluid particles, defined as D{r) oc Ar^, 
here satisfies the famous Richardson scaling law D{r) ^ r^/^jj]. 
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In the diagrams of Figs. [8]d and [8t, skewness and flatness of Am^ are shown in terms 
of r for i = and 0.6. The skewness, H3 is first calculated with Eq. fl37j) . then H4 has 
been determined using Eq. fl48p . At t = 0, l-f^al starts from 3/7 at the origin with small 
slope, then decreases until to reach small values. also exhibits small derivatives near the 
origin, where H4 ^ 3, thereafter it decreases more rapidly than {H^l. At i =0.6, the diagram 
importantly changes and exhibits different shapes. The Taylor scale and the corresponding 
Reynolds number are both diminished, so that the variations of H3 and H4 are associated 
to smaller distances, whereas the flatness at the origin is slightly less than that at t = 0. 
Nevertheless, these variations correspond to higher r/Ay than those for t = 0, and also in 
this case, H4 reaches the value of 3 more rapidly than tends to zero. 

The PDFs of Aur are calculated with Eqs. ([S2D and (jUj), and are shown in Fig. [H] in 
terms of the dimensionless abscissa 

AUr 
^ " (Am2)1/2 

where, these distribution functions are normalized, in order that their standard deviations 
are equal to the unity. The figure represents the distribution functions of s for several r/Ar, 
at t = 0, 0.5 and 0.6, where the dotted curves represent the gaussian distribution functions. 
The calculation of H^{r) is first carried out with Eq. (1571) . then the function ilj{r,Rx) is 
identified through Eq. fl50l) . and finally the PDF is obtained with Eq. fl52l) . For t = (see 
Fig. [9^) and according to the evolutions of and H4, the PDFs calculated at t/At = 
and 1, are quite similar each other, whereas for r/Ay = 5, the PDF is an almost gaussian 
function. Toward the end of the simulation, (see Fig. I9b and c), the two PDFs calculated 
at r/A-r = and 1, exhibit more sizable differences, whereas for r/A-r = 5, the PDF differs 
very much from a gaussian PDF. This is in line with the plots of H^{r) and Hi{r) of Fig. 

El 

Next, the spatial structure of Au^, given by Eq. (H71) . is analyzed using the previous 
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29|, Aur behaves quite 



results of the simulation. According to the various works 
similarly to a multifractal system, where Aur obeys to a law of the kind Aur{r) ~ r'^ where 
the exponent g is a fluctuating function of space. This implies that the statistical moments 
of Aur{r) are expressed through different scaling exponents C{P) whose values depend on 
the moment order P, i.e. 

(Awf (r)) = Ar«^) (57) 
23 
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FIG. 10: Scaling exponents of longitudinal velocity difference versus the order moment at different 



times. Continuous 
mogorov K41 data_ 
She-Leveque data 



ines with solid symbols are for the present data. Dashed lines are for Kol- 
5]. Dashdotted lines are for Kolmogorov K62 data 27|. Dotted lines are for 



These scaling exponents are here identified through a best fitting procedure, in the interval 
2i < r < At, where the statistical moments of Aur{r) are calculated with Eqs. ( |48|) . Figure 
[To] shows the comparison between the scaling exponents here obtained (continuous lines with 
solid symbols) and those of the Kolmogorov theories K41 (dashed lines) and K62 27 1 
(dashdotted lines), and those given by She-Leveque [28| (dotted curves). At t = 0, the slope 
of C(-P) is about constant, whereas the values of ({P) are very different from those calculated 
by the various authors. This means that, for the chosen initial correlation function, Aur{r) 
behaves like a simple fractal system, where C(-P) P- Again, this result depends on the 
fact that, at the initial times, the energy spectrum is not developed. As the time increases, 
the correlation function changes causing variations in the statistical moments of /S.Ur{r). As 
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result, C{P) gradually diminish and exhibit a variable slope which depends on the moment 
order P, until to reach the situation of Fig. [TOb . where the simulation is just ended. The 
correlation function and the dimensionless moments of Aur{r) are changed, thus the plot of 
C{P) shows that near the origin, C(-P) — -P/3, whereas elsewhere the values of C{P) ^-^'s in 
agreement with the She-Leveque results, confirming that Aur{r) behaves like a multifractal 
system. 

Other simulations with different initial correlation functions and Reynolds numbers have 
been performed, and all of them lead to analogous results, in the sense that, at the end of the 
simulations, the diverse quantities such as Q{r), C and C(-P) are quite similar to those just 
calculated. For what concerns the effect of the Reynolds number, its increment determines 
a wider Kolmogorov inertial range and a smaller dissipation energy rate in accordance to 
Eq. fl66l) . whereas the shapes of the various energy spectrums remain qualitatively unaltered 
with respect to Fig. [61 

In order to study the evolution of the intermittancy vs. the Reynolds number. Table 
1 gives the first ten statistical moments of F{dur/dr). These are calculated with Eqs. 
( 148|) and ( 149|) . for Rx = 10.12, 100 and 1000, and are shown in comparison with those of 
a gaussian distribution function. It is apparent that a constant nonzero skewness of the 
longitudinal velocity derivative, causes an intermittancy which rises with Rx (see Eq. (H71) ). 



Moment 


Rx ~ 10 


Rx = 10^ 


Rx = 10^ 


Gaussian 


Order 


P. R. 


P. R. 


P. R. 


Moment 


3 


-.428571 


-.428571 


-.428571 





4 


3.96973 


7.69530 


8.95525 


3 


5 


-7.21043 


-11.7922 


-12.7656 





6 


42.4092 


173.992 


228.486 


15 


7 


-170.850 


-551.972 


-667.237 





8 


1035.22 


7968.33 


11648.2 


105 


9 


-6329.64 


-41477.9 


-56151.4 





10 


45632.5 


617583. 


997938. 


945 



TABLE I: Dimensionless statistical moments of F(dur/dr) at different Taylor scale Reynolds num- 
bers. P.R. as for "present results". 
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FIG. 11: Dimensionless moments -^4(0) and i^6(0) plotted vs. R\. Continuous lines are for the 
present results. The dashed line is the tangent to the curve of i^4(0) in R\ ~ 10. 
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FIG. 12: Flatness H^{^) vs. Rx- These data are from Ref.[21]. 



More specifically, Fig. [TT] shows the variations of i?4(0) and -^6(0) (continuous lines) in 
terms of i?A, calculated with Eqs. ( I48l) and ( 149|) . with -^3(0) = —3/7. These moments are 
rising functions of R\ for 10 < -Ra ^ 700, whereas for higher R\ these tend to the saturation 
and such behavior also happens for the other absolute moments. According to Eq. ( HHI) . in 
the interval 10 < -Ra ^ 70, Hi^ and ifg result to be about proportional to and -Ra^^i 
respectively, and the intermittancy increases with the Reynolds number until to -Ra ~ 700, 
where it ceases to rise so quickly. This behavior, represented by the continuous lines, 
depends on the fact that i\} ~ v'^, and results to be in very good agreement with the data 
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FIG. 13: Skewness S = -^3(0), Flatness F = -^4(0) and hyperflatness Hq{0) vs. Rx. These data 



are from Ref. 



of Pullin and Saffman 



SOj, for 10 < R\ 5, 100. Figure [TT] can be compared with the data 



collected by Sreenivasan and Antonia 



2l| , which are here reported into Fig. [121 These latter 



are referred to several measurements and simulations obtained in different situations which 
can be very far from the isotropy and homogeneity conditions. Nevertheless a comparison 
between the present results and those of Ref. 2l|] is an opportunity to state if the two 
data exhibit elements in common. According to Ref. 2l|], the flatness monotonically rises 
with Rx with a rising rate which agrees with Eq. fH9|) for 10 < i?A ^ 60 (dashed line, Fig. 
TTl) . whereas the skewness seems to exhibit minor variations. Thereafter, H4 continues to 
rise with about the same rate, without the saturation observed in Fig. [TTl The weaker 
intermittancy calculated with the present theory arise from the isotropy which makes the 
velocity fluctuation a gaussian random variable, while, as seen in sec. IVIl without the 
isotropy condition, the flatness of velocity and of velocity difference can be much greater 
than that of the isotropic case. 



Again, the obtained results are compared with the data of Tabeling et al |22|, |23[], where, 
in an experiment using low temperature helium gas between two counter-rotating cylinders 
(closed cell), the authors measure the PDF of du^jdr and its moments. Also in this case the 
flow can be quite far from to the isotropy condition. In fact, these experiments pertain wall- 
bounded flows, where the walls could importantly influence the fluid velocity in proximity of 
the probe. The authors found that the higher moments than the third order, first increase 
with Rx until to Rx ~ 700, then exhibit a lightly non-monotonic evolution with respect to 
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FIG. 14: Log linear plot of the PDF of du^jdr for different R\. (a): dotted, dashdotted and 
continuous lines are for R\ = 15, 30 and 60, respectively, (b) and (c) PDFs for R\ = 255, 416, 
514, 1035 and 1553. (c) represents an enlarged part of the diagram (b) 

i?A, and finally cease their variations denoting a transition behavior (See Fig. [T3!) . As far as 
the skewness is concerned, the authors observe small percentage variations. Although the 
isotropy does not describe the non-monotonic evolution near Rx = 700, the results obtained 
with Eq. fHTl) can be considered comparable with those of Refs. 



22 



23j, resulting also 



in this case, that the proposed theory gives a weaker intermittancy with respect to Refs. 



22 



23|. 



The normalized PDFs of dur/dr are calculated with Eqs. fl52|) and fHTl) . and are shown 

in Fig. [H] in terms of the variable s, which is defined as 

dur/dr 
{{dur/drff'^ 

Figure [T^ shows the diagrams for Rx = 15, 30 and 60, where the PDFs vary in such a way 
that if3(0) = -3/7. 
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FIG. 15: PDF of dur/dr for = 255, 416, 514, 1035 and 1553. These data are from Ref. [23| 
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FIG. 16: Plot of the integrand s'^F{s) for different R\. Dotted, dashdotted and continuous lines 
are for Rx = 15, 30 and 60, respectively. 



As well as in Ref. |23|, Figs. 4b and 4c give the PDF for Rx = 255, 416, 514, 1035 and 1553, 
where these last Reynolds numbers are calculated through the Kolmogorov function given 
in Ref. [23|], with i^3(0) = —3/7. In particular. Fig. [T4b represents the enlarged region of 
Fig. [IDd, where the tails of PDF are shown for 5 < s < 8. According to Eq. fH7|) . the tails 



ligher Rx, smaller variations 
Fig. [Ht shows that the 



of the PDF rise in the interval 10 < -Ra ^ 700, whereas at 
occur. Although the non-monotonic trend observed in Ref. [2 
values of the PDFs calculated with the proposed theory, for 5 < s < 8, exhibit the same 
order of magnitude of those obtained by Tabeling et al 23|] which are here shown in Fig. [T5l 



Asymmetry and intermittency of the distribution functions are also represented through 
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the integrand function of the 4*'* order moment of PDF, which is Ja{s) = s^-F(s) This 
function is shown in terms of s, in Fig. [161 for Rx = 15, 30 and 60. 

VIII. CONCLUSIONS 

The proposed theory is based on the Landau conjecture which states that the turbulence 
is caused by the bifurcations of the velocity field. 

The obtained results confirm the capability of the proposed theory to describe quite well 
the general properties of the turbulence. These results are here summarized: 

1. The analysis of the bifurcations gives the connection between number of bifurcations, 
length scales and Reynolds number at the onset of the turbulence and allows to deter- 
mine the minimum Taylor-scale Reynolds number for isotropic turbulence. This last 
one is about 10, and, below this value, the isotropic turbulence is not allowed. 

2. The momentum equations written using the referential description allow the velocity 
fluctuation to be expressed by means of the Lyapunov analysis of the kinematics of 
fluid deformation. 

3. The Lyapunov analysis of the relative kinematics equations provides an explanation 
of the physical mechanism of the energy cascade in turbulence. The non-ergodicity of 
d/dt{u ■ u'), due to the fluid incompressibility, make possible that the inertia forces 
transfer the kinetic energy between the length scales without changing the total kinetic 
energy. This implies that the skewness of the longitudinal velocity derivative is a 
constant of the present theory and that the energy cascade mechanism does not depend 
on the Reynolds number. 

4. The Fourier analysis of the velocity difference provides the statistics of Aur- This is a 
non-Gaussian statistics, where the constant skewness of dUf/dr implies that the other 
higher absolute moments increase with the Taylor-scale Reynolds number. 

5. The developed energy spectrums, calculated with the proposed closure of the von 
Karman-Howarth equation, agrees quite well with the Kolmogorov law k~^/^ in a 
given interval of k, which defines the inertial subrange of Kolmogorov. 
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6. For developed energy spectrums, the Kolmogorov function is about constant in a 
wide range of separation distances and its maximum is quite close to 4/5, whereas 
the Kolmogorov constant is about equal to 1.93. As the consequence, the maximal 
finite scale Lyapunov exponent and the diffusivity coefficient, vary according to the 
Richardson law when the separation distance is of the order of the Taylor scale. 

7. The proposed theory also describes very well the multifractality of the velocity dif- 
ference, in the sense that, for developed energy spectrum, the scaling exponents of 
the longitudinal velocity difference, when expressed in terms of the moments order, 
exhibit the characteric shape observed by the various authors. 
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X. APPENDIX 



The von Karman-Howarth equation gives the evolution in time of the longitudinal corre- 
lation function for isotropic turbulence. The correlation function of the velocity components 
is the symmetrical second order tensor Rijij) = iuiu'^^, where Ui and u'^ are the velocity 
components at x and x + r, respectively, being r the separation vector. The equations for 
Rij are obtained by the Navier-Stokes equations written in the two points x and x + r p, l7|. 
For isotropic turbulence Rij can be expressed as 



Rijir) = 



{f-9)'-^ + 9h (58) 



/ and g are, respectively, longitudinal and lateral correlation functions, which are 

^ K(x)M^(x + r)) ^^^^ ^ (wn(r)Mn(x + r)) 

where Ur and m„ are, respectively, the velocity components parallel and normal to r, whereas 
r = |r| and = {uf) ={uD= 1/3 {uiUi). Due to the continuity equation, / and g are linked 
each other by the relationship 

S = f^l%r (60) 
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The von Karman-Howarth equation reads as follows 



aa 



K + 2uu' 



+ 



dt \ dr^ r dr ^ 

where K is an even function of r, which is defined by the following equation 

+ 3) K(r) = A {uy,(u, - <)) 



and which can also be expressed as 



d 



where k is the longitudinal triple correlation function 

(m,2(x)u^(x + r)) 



k{r) 



m 



The boundary conditions of Eq. fl6T]) are 

/(0) = 1, lim/(r) = 



(61) 



(62) 



(63) 



(64) 



(65) 



The viscosity is responsible for the decay of the turbulent kinetic energy, the rate of which 
is obtained putting r = in the von Karman-Howarth equation, i.e. 



dv? 



(0) 



(66) 



dt dr"^ 

This energy is distributed at different wave-lengths according to the energy spectrum E{k,) 
which is calculated as the Fourier Transform of fu^, whereas the "transfer function" T(k) 

n 

is the Fourier Transform of LO, i.e. 



EiK) 








1 ("^ 




T^Jo 




T{k) 




K{r) 



2 2 

K r 



sm K,r 



nr 



cos nr 1 dr 



(67) 



where k = \n\ and T(k) identically satisfies to the integral condition 

/■oo 

/ T{K)dK = 
Jo 



(68) 



which states that K does not modify the total kinetic energy. The rate of energy dissipation 
e is calculated for isotropic turbulence as follows 

3du^ 



2 dt 



POO 

2V I K^E{K)dK 

Jo 



(69) 
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The microscales of Taylor At, and of Kolmogorov I, are defined as 
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